Thermodynamically admissible form for discrete hydrodynamics 
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We construct a discrete model of fluid particles accord- 
ing to the GENERIC formalism. The model has the form of 
Smoothed Particle Hydrodynamics including correct thermal 
fluctuations. A slight variation of the model reproduces the 
Dissipative Particle Dynamics model with any desired ther- 
modynamic behavior. The resulting algorithm has the follow- 
ing properties: mass, momentum and energy are conserved, 
entropy is a non-decreasing function of time and the thermal 
fluctuations produce the correct Einstein distribution func- 
tion at equilibrium. 

Particle based methods for solving hydrodynamic 
problems are good candidates for the study of complex 
fluids because they allow an easy treatment of compli- 
cated geometries as those appearing in the interstitial re- 
gions of a colloidal or polymeric suspension. In addition, 
they allow us to use molecular dynamic codes which are 
comparatively much simpler than usual computational 
fluid dynamics algorithms. 

Two seemingly distinct particle methods are partic- 
ularly appealing because of their versatility, Smoothed 
Particle Hydrodynamics (SPH) and Dissipative Particle 
Dynamics (DPD). SPH is a well-known model in the con- 
text of computational astrophysics jjj that is receiving 
growing interest for the study of laboratory fluid dynam- 
ical problems ||. This method is basically a discretiza- 
tion of Navier-Stokes equations on a Lagrangian grid with 
the aid of a weight function. At present, there is still no 
version of SPH that consistently includes thermal fluc- 
tuations, that is, there is no SPH discretization of fluc- 
tuating hydrodynamics However, thermal fluctua- 
tions are crucial if one wants to use SPH for the study 
of complex fluids. Thermal fluctuations are the ultimate 
responsible for the diffusive behavior of suspended ob- 
jects ||. We are faced, therefore, with the problem of 
generalizing SPH to the mesoscopic level where fluctu- 
ations are important. On the other hand, DPD || is 
another particle based model that includes thermal fluc- 
tuations in a sensible way [Tj and it has hydrodynamic 
behavior ||. For these reasons it has been applied to 
a large variety of problems dealing with complex fluids 
The initial problem of non-conservation of the en- 
ergy has been resolved and the technique can now be 
applied to non-isothermal problems as well JlO| . There 
remains, however, a fundamental problem in DPD which 
is the physical interpretation of the conservative interac- 
tions between DPD particles, which determine the full 
thermodynamic behavior of the system |Tl]] . 



In this letter, we give a solution to these basic prob- 
lems with SPH and DPD. We formulate a Lagrangian 
discrete off-lattice model for hydrodynamics which gen- 
eralizes SPH by including correct thermal fluctuations. 
A slight variation produces the DPD model with any de- 
sired thermodynamic behavior. The derivation of these 
models is done within the context of the GENERIC for- 
malism which ensures thermodynamic consistency [ p"2| . 
GENERIC postulates that any physically sensible dy- 
namic equation in non-equilibrium thermodynamics has 
an underlying structure. All the dynamic equations stud- 
ied so far fit into the formalism. Linear irreversible ther- 
modynamics, non-relativistic and relativistic hydrody- 
namics, Boltzmann's equation, polymer kinetic theory, 
and chemical reactions, just to mention a few, have all the 
GENERIC structure It is natural, then, to for- 

mulate a model for fluid particles within the GENERIC 
formalism. As it is becoming apparent in recent years, 
by putting "more physics" into discretization of partial 
differential equations one gets improved behaviour (i.e. 
stability) of the discretized equations. 

The GENERIC formalism states that the time evolu- 
tion equations for a complete set of independent variables 
x required for the description of a non-equilibrium system 
have the structure 
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The first term in the right hand side produces the re- 
versible part of the dynamics whereas the second term is 
responsible of the irreversible dissipative dynamics. Here, 
E, S are the energy and entropy of the system expressed 
in terms of the variables x and L,M are matrices that 
satisfy the following degeneracy requirements, 
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dE 
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dx 
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In addition, L is antisymmetric (this guarantees that en- 
ergy is conserved) and M is a positive definite symmetric 
matrix (this guarantees that the entropy is a nondecreas- 
ing function of time). If the system presents dynamical 
invariants I(x) different from the total energy, then fur- 
ther restrictions on the form of L and M are required 
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These conditions ensure that dl/dt = 0. 
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The deterministic equations ([[]) are, actually, an ap- 
proximation in which thermal fluctuations are neglected. 
If thermal fluctuations are not neglected, the dynamics is 
described by stochastic differential equations or, equiv- 
alently, by a Fokker-Planck equation that governs the 
probability distribution function p = p(x,t). This FPE 
has the form E^] 



d t p = - 

Ox 
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where k B is Boltzmann's constant. We require that the 
equilibrium solution of this Fokker-Planck equation is 
given by the Einstein distribution function in the pres- 
ence of dynamical invariants [|14|, this is 



p°*{x) = g(E(x), I(x)) e X p{S(x)/k B }, 



(5) 



where the function g is completely determined by the ar- 
bitrary initial distribution of dynamical invariants. This 
imposes further conditions on the matrices L, M, namely, 
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(6) 



The first property can be derived independently with pro- 
jection operator techniques fj^| . The second property 
implies that the last equation in (j^) is automatically sat- 
isfied. When fluctuations are present, the entropy func- 
tional S[p] — J S(x)p(x,t)dx — k B J p{x, t) In p(x, t)dx 
plays the role of a Lyapunov function with > 0. 

We finally consider the Ito stochastic differential equa- 
tions that are mathematically equivalent to the above 
Fokker-Planck equation ]15[ 



dx 
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dt + dx. 
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Here, the stochastic term dx is a linear combination of in- 
dependent increments of the Wiener process. It satisfies 
the mncmotcchnical Ito rule 



dxdx T = 2k B Mdt, 



(8) 



which means that dx is an infinitesimal of order 1/2 
jL5|. Eqn. (|J) is a compact and formal statement of 
the fluctuation-dissipation theorem. 

When formulating new models it might be convenient 
to specify dx directly instead of M. This ensures that 
M through (^) automatically satisfies the symmetry and 
positive definite character. In order to guarantee that 
the total energy and dynamical invariants do not change 
in time, a strong requirement on the form of dx holds, 



dE , dl , 

-7T— -dx = 0, —-dx — 0, 

Ox Ox 

implying the last equations in (||) and 



(9) 



The basic problem in any non-equilibrium description 
of a given system is to identify the relevant set of variables 
in the system. We aim at modeling fluid particles, that is, 
portions of fluid or large clusters of molecules that move 
following collective motions. It is sensible to assume that 
these fluid particles are small moving thermodynamic 
systems with the center of mass located at r,, and with 
momentum p^, volume Vi, mass rrii, and entropy Sj (or in- 
ternal energy e,). Even though they are thermodynamic 
systems in the sense that an entropy function can be de- 
fined, the fluid particles are assumed to be small enough 
to suffer from stochastic fluctuations due to the underly- 
ing molecules forming the fluid particle. The state of the 
system is x = {r^, p^, Vi, m^, s^, i = 1, . . . , N}, where N 
is the number of fluid particles. 

The two basic building blocks in GENERIC are the 
energy and the entropy of the system as a function of the 
selected variables. In our case they are 



where e(Vi, nij, Si) is the internal energy as a function of 
the extensive variables of the fluid particle. We will as- 
sume the hypothesis of local equilibrium, in accordance 
with the usual treatment of hydrodynamics. Therefore, 
the internal energy of the fluid particles is the same func- 
tion of the mass, volume, and entropy as the total internal 
energy of the whole fluid system at equilibrium. We can 
now consider the derivatives of E and S which are given 
by 
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(ii) 



where = pi/rrn is the velocity of the fluid particles and 
the intensive parameters are the pressure pi, the chem- 
ical potential per unit mass pi, and the temperature Tj 
which arc functions of the extensive variables Vi,rrii, Si of 
the fluid particles. We have to consider also the dynami- 
cal invariants I(x) that we want to retain in the discrete 
model. These dynamical invariants are the total momen- 
tum P = YliPi> tne total volume V = J^Vj, and the 
total mass M = ^ m i- 

We construct now the L matrix for the discrete hy- 
drodynamics problem. We impose that the reversible 
part of the equations of motion should give r.i| rev = Vj, 
?Tii|rev = Sj|rcv = 0. That is, the position of the fluid 
particles changes according to its velocity, the "identity" 
of the fluid particles is given by the constant mass it pos- 
sesses, and finally, the reversible part of the dynamics 
should not produce any change in the entropy content of 
the fluid particle. By looking at the term L-dE/dx in 
Eqn. ^ with the energy derivatives given in ( pi] ) , the L 
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matrix that produces the desired equations can only be 
of the form of N x N blocks Ly of size 9 x 9 of the form 



(12) 



The two last rows ensure the invariance of rrii,Si, the 
two last columns are fixed by antisymmetry. The first 
row ensures the equation of motion for the position and 
the first column is fixed by antisymmetry. We still have 
freedom for the form of the vector fly which can, in prin- 
ciple, depend on all state variables. L is antisymmetric 
because Ly = — and it satisfies Eqn. (0). The effect 
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on the invariants, Eqn. 



|3j), is guaranteed if Hy obeys 

E n ii=°' ( l3 ) 



Finally, we can guarantee the first property in Eqn. 
if the following identity is satisfied 



£ 



d Pi ijPj + dVi 



0. 



(14) 



The resulting reversible equations for momentum and 
volume are 



Vi 



(15) 



Note that the vector fiy can be interpreted as a discrete 
version of the gradient operator in such a way that the 
above equations look like a discrete version of the mo- 
mentum equation and continuity equation (in terms of 
the specific volume 1/p instead of the mass density p) of 
a non-dissipative fluid in a Lagrangian description. We 
propose the following form for fly 



1 

N 



E^-E 



(16) 



where ajy = (V/A^) 2 VA(ry), and A(r) being a weight 
function of range r c normalized to unity, J drA(r) = 1. 
This form @ satisfies Eqns. @ and @. If the range 
r c is much smaller than the typical length scale of vari- 
ation of the hydrodynamic variables and the typical dis- 
tance between points is much smaller than r c , then it is 
possible to prove that 



v<v/(r,-)~x>«/fa) 



(17) 



Under these circumstances, the discrete equations ( |15| ) 
converge towards the continuum Euler equations of hy- 
drodynamics. 



Our aim now is to construct the matrix M. We have 
to specify first where irreversibility occurs. We do not 
want irreversible processes associated with the evolution 
of the position, volume or mass of the particles. This 
implies that the noise term in the equation of motion (j?]) 
has the structure dx T — » (0, dpi, 0, 0, dSi). 

Thermal fluctuations are introduced into the equations 
of hydrodynamics through a random stress tensor and 
random heat flux [0. By simple analogy with the ex- 
pression of the random noise in non-linear hydrodynam- 
ics pi| , we postulate the following random terms 

dpi = y^Sljj-do-j, 

3 



3 3 

where the random stress d&i and random heat flux d3\ 
are defined by 

da t = {4k B T im ) 1 / 2 ~dW- + (ekBTibf'HjjixldWi], 

(19) 



d~3l = T l (2fc sKl ) 1/2 dV, 



Here, r\i is the shear viscosity, Q is the bulk viscosity, and 
Ki is the thermal conductivity. The traceless symmetric 

random matrix cfW \ is given by 

5Wf = \ [dWi + dWf] - ltr[dW,]l. (20) 

D is the physical dimension of space and dWi is a matrix 
of independent Wiener increments. The vector dV \ is 
also a vector of independent Wiener increments. They 
satisfy the Ito mnemotechnical rules 



dWf^ dW u , 



0. 



(21) 



Note that the postulated forms for dpi, dSi in Eqn. 
satisfy J2i v «'rfPi + Tidsi — and J2i dpi — and, there- 
fore, Eqns. (^) are satisfied. This means that the pos- 
tulated noise terms conserve momentum and energy ex- 
actly. It is now a matter of algebra to construct the 
dyadic dxdx T and from Eqn. (||) extract the matrix M. 
Once M is constructed, the final equations of motion for 
the discrete hydrodynamic variables are (D = 3 and, for 
simplicity, constant transport coefficients are assumed), 

drrii = 0, dti — v;di, dVi = Didt 

dpi = £ ^iiPi + E dt + ^ 



C v 

^Slij-J^dt + Tids. 



2r)Gi : Gi + £D*) dt 



(22) 
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We have introduced the stress tensor <x,-, the traceless 
symmetric velocity "gradient" tensor Gj, and the "di- 
vergence" Di by 

of = -2r?Gr ~ CA<^ 

of = \ Y\^i + n^ k ] ]T -v* 

fe ft 

A = ^ ^ifc-Vfe 



The heat flux of particle i is defined by, 



J? 



1 



Cvk 



(l + 5 lk ) 



(23) 



(24) 



One can easily recognize all the terms corresponding to 
the continuum equations of hydrodynamics ||l7j] . The fac- 
tor ks/Cvi, where Cv% is the specific heat at constant 
volume of particle i, comes from the term ksd-M/dx in 
Eqn. (0). In the continuum version of hydrodynamics, 
this term is zero due to locality |jq| . 

Eqns. (p^ ) are the main result of this Letter. They 
have the structure of SPH but conserve total mass, total 
momentum, total energy and total volume of the parti- 
cles. Because of the GENERIC structure of the equa- 
tions, the entropy functional S[p] of the system is a non- 
decreasing function of time and the equilibrium solution 
is given by Einstein distribution function. 

Instead of the noise terms dlq) one could also postulate 
the noise terms as in DPD ]10[ 

dpi = ^BjjdWjj 



dsi 



y E B • V <' /1, '<.- • fx A v dV « ( 25 ) 



where Bjj = — Bjj and = Aj$ are suitable func- 
tions of position and, perhaps, other state variables. The 
independent Wiener processes satisfy dWij — dWji and 
dVij = —dVji and the following Ito mnemotechnical rules 

dWa'dWjj' = [SijSi'j' +5 l j'8i'j}dt 
dVu'dVjj, = [SijS^f - SifSi/j]dt (26) 

These noise terms satisfy also Eqn. (^) and momentum 
and energy are conserved. The resulting dynamic equa- 
tions have the same reversible part as in Eqn. (p2|) and 
have the familiar "Brownian dashpot" dissipative forces 
of the DPD model. In this way, by introducing a vol- 
ume and an internal energy variables into the standard 
DPD model, one derives a thermodynamically consis- 
tent model in which the "conservative" forces are truly 
pressure forces between DPD particles. The resulting 
DPD equations are simpler than Eqns. ( p2| ) and produce 
macroscopic hydrodynamic behaviour, but the identifi- 
cation of the actual values of the transport coefficients is 
less obvious. 
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